function fai(t::Float64,t0::Float64)
    if t < t0
        println("error: t should be greater than t0")
        return nothing
    end
    RH = 0.5
    h = 184.0
    fcm = 38
    beta_H = 150.0*(1 + 1.2*(RH)^18)*h/100.0 + 250.0
    if beta_H > 1500.0
        beta_H = 1500.0
    end
    beta_ctt0 = ((t - t0)/(beta_H + t - t0))^0.3
    v1 = 1 + (1 - RH)/(0.46*(h/100.0)^(1/3))
    v2 = 5.3/(fcm/10.0)^0.5
    v3 = 1.0/(0.1 + t0^0.2)
    phi_0 = v1 * v2 * v3

    return phi_0*beta_ctt0
end

function j_fun(t::Float64,t0::Float64)
    Eci = 34000.0
    s = 0.25
    beta_cct = exp(s*(1.0 - sqrt(28.0/t0)))
    Ecit0 = sqrt(beta_cct) * Eci
    return 1.0/Ecit0 + fai(t, t0)/Eci
end

function epsilon_sh(t::Float64, ts::Float64)
    RH = 0.5
    beta_sc = 5.0
    fcm = 38.0
    h = 184.0
    beta_SRH = 1.0 - RH^3;
    beta_RH = 0.4 <= RH < 0.99 ? -1.55*beta_SRH : 0.25
    epsilon_fcm = (160 + 10*beta_sc*(9.0 - fcm/10.0))*1e-6
    epsilon_cs0 = epsilon_fcm*beta_RH
    beta_stt0 = sqrt((t - ts)/(350*(h/100.0)^2 + t - ts))
    return epsilon_cs0*beta_stt0
end

function beta_E(t::Float64)
    s = 0.25
    return sqrt(exp(s*(1.0 - sqrt(28.0/t))))
end

function j_matrix(i::Int64, j::Int64, time::Vector{Float64})
    return j_fun(time[i], time[j])
end

function en(n::Int64, time::Vector{Float64})
   if n == 1
       r = 1.0/j_matrix(1, 1, time)
   elseif n >= 2
       r = 2.0/(j_matrix(n, n, time) + j_matrix(n, n - 1, time))
   else
       println("error: n is illeagal")
       return nothing
   end
   return r
end

function psi(n::Int64, i::Int64, time::Vector{Float64})
    if i == 1
        r = (j_matrix(n, 2, time) - j_matrix(n, 1, time))/(j_matrix(n, n, time) + j_matrix(n, n - 1, time))
    elseif i > 1
        r = (j_matrix(n, i + 1, time) - j_matrix(n, i - 1, time))/(j_matrix(n, n, time) + j_matrix(n, n - 1, time))
    end
    return r
end

function epsilon_shn(n::Int64, time::Vector{Float64})
    IsShrink = 1.0
    return IsShrink*epsilon_sh(time[n], time[1])
end


#=
time0 = [7.0, 30.0, 50.0, 80.0, 120.0, 170.0, 230.0, 300.0, 365.0]
fai(365.0, 7.0)
j_fun(365.0, 7.0)
epsilon_sh(150.0,7.0)
en(7,time0)
psi(8,3,time0)
j_matrix(4,3,time0)
=#
